library(HEIfunctions)
library(spBayes)
library(splines)
library(xtable)
library(pscl)
library(DMwR)
set.seed(51)

load('datforimpute_withtemp.Rda')

dat$baseorig1990_1994 = dat$pm10base1990_1994
dat$baseorig1990 = dat$pm10base1990

#################################################################################
##########    Merge Baseline Predictions for pm10base1990   ##############
#################################################################################
load(file = "20150427pm1990preds.Rda")
####        Check Impute Baseline Predictions.R for nsamp and nburn.  Burnin has already been taken out to calculate outmodel$pred, but not for the raw samples (e.g., outmod$model$p.theta.samples)            ######
nsamp    <- 50000
nburn		<- 10000
obsdat = subset(dat, !is.na(pm10base1990))
misdat = subset(dat, is.na(pm10base1990))

#Grab the predicted mean, sd, and 5 individual samples from the posterior-predictive distribution (for possible sensitivity analysis later)
whichpostpreds = sample(nburn:nsamp, size=5)
imputations = cbind(outmodel$pred[, -1], outmodel$predmodel$p.y.predictive.samples[, whichpostpreds])
impdat = as.data.frame(cbind(as.character(outmodel$pred$Monitor), imputations))
names(impdat) = c("Monitor", "predmean1990", "predsd1990", paste("ppbase1990", 1:length(whichpostpreds), sep=""))

datmerged = merge(dat, impdat, by="Monitor", all.x=TRUE)
with(datmerged, table(is.na(baseorig1990), is.na(predmean1990)))
datmerged$pm10base1990[is.na(datmerged$pm10base1990)] = datmerged$predmean1990[is.na(datmerged$pm10base1990)]

dat = datmerged

### Output 
save(dat, file = 'pm10dat_withbaselineimp.RData')





